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Abstract 

The alternating-current optimal power flow (ACOPF) is one of the 
best known non-convex non-linear optimisation problems. We present a 
novel re-formulation of ACOPF, which is based on lifting the rectangu¬ 
lar power-voltage rank-constrained formulation, and makes it possible to 
derive alternative SDP relaxations. For those, we develop a first-order 
method based on the parallel coordinate descent with a novel closed-form 
step based on roots of cubic polynomials. 

1 Introduction 

Alternating-current optimal power flow problem (ACOPF) is one of the best 
known non-linear optimisation problems [59) . Due to its non-convexity, de¬ 
ciding feasibility is NP-Hard even for a tree network with fixed voltages [50] . 
Still, there has been much recent progress [551 [50]: Bai et al. [5] introduced 
a semidefinite programming (SDP) relaxation, which turned out to be partic¬ 
ularly opportune. Lavaei and Low [29] have shown that the SDP relaxation 
produces exact solutions, under certain spectral conditions. More generally, it 
can be strengthened so as to obtain a hierarchy of SDP relaxations whose op¬ 
tima are asymptotically converging to the true optimum of ACOPF [l9l [251 [22] . 
Unfortunately, the run-times of even the best-performing solvers for the SDP 
relaxations [551 EH remain much higher than that of commonly used methods 
without global convergence guarantees such as Matpower [5^ ■ 

Following a brief overview of our notation, we introduce a novel rank-constrained 
reformulation of the ACOPF problem in SectionEl where all constraints, except 
for the rank constraint, are coordinate-wise. Based on this reformulation, we 
derive novel SDP relaxations. Next, we present a parallel coordinate-descent 
algorithm in Section |4l which solves a sequence of convex relaxations with 
coordinate-wise constraints, using a novel closed-form step considering roots 
of cubic polynomials. We can show: 


1 


• the algorithm converges to the exact optimum of the non-convex problem, 
when the optimum of the non-convex problem coincides with the opti¬ 
mum of the novel SDP relaxations and certain additional assumptions are 
satisfied, as detailed in Section [5] 

• the algorithm suggests a strengthening of the novel SDP relaxations is 
needed, whenever it detects the optimum of the non-convex problem has 
not been found, using certain novel sufficient conditions of optimality 

• our pre-liminary implementation reaches a precision comparable to the de¬ 
fault settings of Matpower [62) , the commonly used interior-point method 
without global convergence guarantees, on certain well-known instances 
including the IEEE 118 bus test system, within comparable times, as de¬ 
tailed in Section [3 

although much more work remains to be done, especially with focus on large- 
scale instances. Also, as with most solvers, the proposed assumes feasibility and 
does not certify infeasibility, when encountered. The proofs of convergence 
rely on the work of Burer and Monteiro [12] , Grippo et al. [20], and our earlier 
work piHT] . 


2 The Problem 

Informally, within the optimal power flow problem, one aims to decide where to 
generate power, such that the demand for power is met and costs of generation 
are minimised. In the alternating-current model, one considers the complex volt¬ 
age, complex current, and complex power, although one may introduce decision 
variables representing only a subset thereof. The constraints are non-convex in 
the alternating-current model, and a particular care hence needs to be taken 
when modelling those, and designing the solvers to match. 

Formally, we start a network of buses N, connected by branches L C N x N, 
modeled as Il-equivalent circuits, with the input comprising also of: 

• G C N, which are the generators, with the associated coefficients 
of the quadratic cost function at /c S G, 

• Pk + jQk: which are the active and reactive loads at each bus k G N, 

• P™™, Q/f™ and which are the limits on active and reactive 

generation capacity at bus k G G, where = 

0 for all A: G iV \ G, 

• yG Cl'^lxl^l, which is the network admittance matrix capturing the value 
of the shunt element bim and series admittance gim + jhm at branch 
(?, m) G L, 

• be™™ and which are the limits on the absolute value of the voltage 

at a given bus k, 
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• which is the limit on the absolute value of the apparent power of a 
branch {l,m) S L. 

In the rectangular power-voltage formulation, the variables are: 

• + jQfc) which is the power generated at bus k G G, 

• Pim + jQim, which is the power flow along branch {I, m) G L, 

• 5RI4 + j^Vk, which is the voltage at each bus k G N. 

The power-flow equations at generator buses k G G are: 

n n 

Pk = Pk + ^Vk + ^Vk + ^y^km), (i) 

2=1 2=1 

n n 

Ql = Qt + ^Vk Yi-’^y^km - ^y^k‘^V,) + ^Vk 

2=1 2=1 

( 2 ) 


while at all other buses k G N\G we have: 

n n 

0^P^ + ^kYi^y^^m - ^y^k^Vi) + 314 ^(32/,fe5?I4 + Jiy.fcSV,), (3) 

2=1 2=1 

n n 

0^Qi+ KFfc - Uy^k^Vi) + - 3y,fc314)- (4) 

2=1 2=1 

Additionally, the power flow at branch (/,m) € L is expressed as 
Plm = bimim^Vm - _ SR^^SR^^), 

(5) 

Qim = bUm^Vm - 31^314. - - 3^,2) 

+ gimrnVi^Vm - 5^K^3y^ - KKn3I^) - (6) 

Considering the above, the alternating-current optimal power flow (ACOPF) is: 

min {cl{Pi + Pif + cl{Pi + Pi) + cl) [ACOPF] 

keG 

s.t. Pf” < Pi + Pi < Pr^ WkGG 

Qt"<Ql + Qt<Qr^ VfceG 

(^mm)2 < sjjy^2 <^y^2 < (ymax)2 'ik G N 

PL + QL<{STrin^ y{l,m)GL 
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where is the coefficient of the leading term of the quadratic cost function at 
generator k. 

In line with recent work [IIMIIIIIIS], let Ck be the standard basis vector 
in and dehne: 


Vk 

yim 

Yk 

Yk 

Mk 

^Im 

^Im 


ekcly, 


(j 




■ gim + jbim)eief - {gim + jbim)eie. 


T 
m 5 


1 

2 


^ivk +yl) 
^ivk - yl) 


S(2/fc - Vk) 

"^{vk + yD. 


1 

2 


s(yfc + Vk) 

^Ivl - Vk) 


^{vk - yl) 
‘^{vk + yl) 


euel 0 
0 etel 


1 

2 

-1 

'Y 


K(y/™ + ylY) ’^{yln - yin.) 
^{yim - yin) ^{yim + yln)\ ’ 
^{yim + yin) ^{yim - yin) 

. ^ivln - yim) ^{yim + vln) 


Using these matrices, we can rewrite [ACOPF] as a real-valued polynomial opti¬ 
mization problem of degree four in variable x € comprising of 5iI4 G 

and 314 S R^^l, stacked: 


™ (4(ir(YkXx^) + -Pfc + c^(tr(Yfea;a;'^) -h P^) + cl) 


keG 


pmax 

k 


s.t. pr < triY.xx^ ) + Pl< PI 
PI = triXkXx'^) 

Qf" < tT{Y,xx^) + Qt< Qr" 
Qi = triYkXx'^) 

< tr(Mfe:r:r^) < 

{tT{YiniXx'^))^ + {tT{YiniXx'^))'^ < 


[PP4] 

VfceG (7) 
VfceA\G (8) 
VfceG (9) 
Vfc e A\G (10) 
VkG N (11) 
V(?, in) G L, (12) 


Henceforth, we use n to denote 2|iV|, i.e., the dimension of the real-valued 
problem |[PP4][ 

Further, the problem can be lifted to obtain a rank-constrained problem in 
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W = xx'^ € 




[Rl] 

keG 



< tr(rfcW) + < Pf 

\/kGG 

(13) 

P^ = triYkW) 

\/kGN\G 

(14) 

Qr’^<tT{YkW) + Qt<Qr^ 

WkGG 

(15) 

Qt = tv{YkW) 

Vfc e iv \ G 

(16) 

(ymi„)2 < tr(MfeW) < ( 14 “'"’')^ 

VkGN 

(17) 

{iv{Yi^W)f + {iv{Yi^W)f < 

V(Z, m) G L 

(18) 

IT P 0, rank(lT) = 1, 


(19) 


for a suitable definition of fk ■ This problem |[R1]| is still NP-Hard, but where 
one can drop the rank constraint to obtain a strong and efficiently solvable SDP 
relaxation: 


min^/fe(lT) s.t. dT^-TTa. WhO, 


[LL] 


keG 


as suggested by [5]. Lavaei and Low [! 
the relaxation [LL] (Optimization 3 of 


studied the conditions, under which 
]) is equivalent to [Rljl We note that 


for traditional solvers the dual of the relaxation (Optimization 4 of 

[2^ 1 is much easier to solve th an|[LL][ as documented in Table [2j Ghaddar et al. 
[19] have shown the relaxation [LL] to be equivalent to a first-level [OP 4 — Hi] of 
a certain hierarchy of relaxations, and how to obtain the solution to |[Rl ] under 
much milder conditions than those of 1291. 


3 The Reformulation 

The first contribution of this paper is a lifted generalisation of |[R1][ 
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min^ fk(W) 


[RrBC] 

keG 



s.t. tk = tr(YfcTD) 

ykGN 

(20) 

pf” - p^ <tk< pr^ - Pi 

WkGG 

(21) 

Pi = tk 

Vfc G iV \ G 

(22) 

gu = tT{YkW) 

WkGN 

(23) 

Qf" -Qt<9k< Qr^ - Qt 

ykGG 

(24) 

Qk =9k 

ykGN\G 

(25) 

hk = tT{MkW) 

WkGN 

(26) 

(t4““)2 <hk< ( 14 ““)" 

WkGN 

(27) 

'^im — tr(TTyj.jTT) 

y{l,m) G L 

(28) 

'^Im — tr(L/mTE) 

y{l,m) G L 

(29) 

— {P'lm) P {vim) 

y{l,m) G L 

(30) 

^ ^ ( cmax\2 

^Irfi ^ Wlm ) 

\/{l,m) G L 

(31) 

W P 0,rank(bE) < r. 


(32) 


There, we still have: 

Proposition 3.1 (Equivalence). [RIBC] is equivalent to\[PP4]\ 


Even in the special case of r = 1, however, we have lifted the problem to a 
higher dimension by adding variables t, g,h,u,v, z, which are box-constrained 
functions of W. 

Subsequently, we make four observations to motivate our approach to solving 


the [RrBCl 


Oi: constraints (HH), (El, (El, and El are box constraints, while the re¬ 
mainder of (I 2 OH 3 T]) are linear equalities 


O 2 ' using elementary linear algebra: 


{W G 5" s.t. W PQ, rank(bE) < r} 
= {RR^ s.t. R G R”"'!, 


where 5" is the set of symmetric n x n matrices. 


O 3 : if rank(lE*) > 1 for the optimum W* of [LL][ there are no known methods 
for extracting the global optimum of [Rl] from W, except for |19) . 


O 4 : zero duality gap at any SDP relaxation in the hierarchy of [T9] does not 

c.f. ETj. 


guarantee the solution of the SDP relaxation is exact for [PP4 


Note that Lavaei and Low [29] restate the condition in Observation O 3 in 
terms of ranks using a related relaxation (Optimization 3), where the rank has 
to be strictly larger than 2 . 
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4 The Algorithm 

Broadly speaking, we use the well-established Augmented Lagrangian approach 
isam], with a low-rank twist |12j . and a parallel coordinate descent with a 
closed-form step. Considering Observation O 2 , we replace variable W G 
by RR^ S p^rtxn for R G to obtain the following augmented Lagrangian 

function: 

£(i?,t,/i,5,w,u,z,A‘,A9,A^,A“,A",A^) := (33) 

E 

keG 

- E ^kitk - tviYkRR^)) + ^ E(^fc - tT{YkRR^)f 

k k 

- E ^'k(9k - triXkRR^)) + f E(5fc - tr(yfeRR^))2 

k k 

- E ^kihk - tr(Mfci?i?^)) + ^Y.^hk- tr{MkRR^)f 

k k 

- E - ty{Y(l^m)RR^)) + f E 

{l,m) 

- E + f E “ H%,m)RR^))‘^ 

{l,m) {1-,'tTT') 

(Z,m) (Z,m) 


Note that constants > 0 pre-multiply regularisers, where TZ can often be 
0 in practice, although not in our analysis, where we require TZ = det(i?^i?), 
which promotes low-rank solutions. 

As suggested in Algorithm [ 1 ] we increase the rank r = 1 , 2 ,... allowed in 
W in an outer loop of the algorithm. In an inner loop, we use coordinate descent 
method to find an approximate minimizer of C{R, t, h, g, u, v, z, A‘, A®, A^, A“, A'“, A^), 
and denote the A:-th iterate R^^'>. Note that variables t,h,g,z have simple box 
constraints, which have to be considered outside of C. 

In particular: 


The Outer Loop 

The outer loop (Lines 1-12) is known as the “low-rank method” [T^]. As sug¬ 
gested by Observation O 3 , in the case of [LL] one may want to perform only 
two iterations r = 1, 2. In the second iteration of the outer loop, one should like 
to test, whether the numerical rank of the iterate Rk in the inner iteration k 
has numerical rank 1. If this is the case, one can conclude the solution obtained 
for r = 1 is exact. This test, sometimes known as “flat extension”, has been 
studied both in terms of numerical implementations and applicability by Burer 
and Choi fTO]. 
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Algorithm Schema 1: A Low-Rank Coordinate Descent Algorithm 
1: for r = 1, 2 , ... do 
2 : choose R € 

3: compute corresponding values of t, h, g, u, v, z 

4: project t, h, g, u, v, z onto the box constraints 

5: for fc = 0, 1, 2,. .. do 

6 : in parallel, minimize C in t, g, h, u, v, z, coordinate-wise 

7: in parallel, minimize C in _R, coordinate-wise 

8 : update Lagrange multipliers A*, A®, A^, A“, A’', A^ 

9: update 

10 : terminate, if criteria are met 

11: end for 

12: end for 


The Inner Loop 

The main computational expense of the proposed algorithm is to find an ap¬ 
proximate minimum of £(i?, t,/i, 5 , it, u, z, A‘, A®, A^, A“, A", A^) with respect to 
i?, t,/i, 5 , 11 , u, z within the inner loop (Lines 6-7). Note that £ as a function 
of i? is - in general - non-convex. The inner loop employs a simple iterative 
optimization strategy, known as the coordinate descent. There, two subsequent 
iterates differ only in a single block of coordinates. In the very common special 
case, used here, we consider single coordinates, in a cyclical fashion. This algo¬ 
rithm has been used widely at least since 1950s. Recent theoretical guarantees 
of random coordinate descent algorithm are due to Nesterov [33] and the present 
authors [Ml ED- See the survey of Wright m for more details. 

The Closed-Form Step 

An important ingredient in the coordinate descent is a novel closed-form step. 
Nevertheless, if we update only one scalar of i? at a time, and fix all other scalars, 
the minimisation problem turns out to be the minimisation of a fourth order 
polynomial. In order to find the minimum of a polynomial ax'^+bx^+cx'^+dx+0, 
we need to find a real root of the polynomial Aax^ -I- Zhx^ -I- 2cx -I- d = 0. This 
polynomial has at most 3 real roots and can be found using closed form formulae 
due to Cardano [13]. Whenever you fix the values across all coordinates except 
one, finding the best possible npdate for the one given coordinate requires either 
the minimisation of a quadratic convex fnnction with respect to simple box 
constraints (for variables t,g,h,z) or minimisation of a polynomial function of 
degree 4 with no constraints (for variables R, u, v), either of which can be done 
by checking the objective value at each out of 2 (for variables t, g, h, z) or 3 (for 
variables R, u, v) stationary points and choosing the best one. 
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The Parallelisation 

For instances large enough, one can easily exploit parallelism. Notice that min¬ 
imization of coordinates of t, g, h, u, v, z can be carried out in parallel without 
any locks, as there are no dependences. One can also update coordinates of R 
in parallel, although some degradation of the speed-up thus obtainable is likely, 
as there can be some dependence in the updates. The degradation is hard to 
bound. Most analyses, c.f. m, hence focus on the uniformly random choice of 
(blocks of) coordinates, although there are exceptions [H]. Trivially, one could 
also parallelise the outer loop, for each r that should be considered. 

Sufficient Conditions for Termination of the Inner Loop 

For both our analysis and in our computational testing, we use a “target infea¬ 
sibility” stopping criterion for the inner loop, considering squared error: 

=Y.{tu - tr(ni?«(i?W)^))2 + Y,{9k - tr(ni?('=)(i?W)^))^ + 

k k 

(34) 

Y,{hk - tr(Mfci?W(i?«)^))2 + ^ («(;_„) - tr(y(i,„)i?W(i?W)^))2+ 

k 

E + E - Am) ” vlum)f ■ 

{l,m) 

We choose the threshold to match the accuracy in terms of squared error ob¬ 
tained by Matpower using default settings on the same instance. Tfe(i?(^)) < 
0.00001 is often sufficient. 

The Initialisation 

In our analysis, we assume that the instance is feasible. This is difficult to 
circumvent, considering Lehmann et al. [30] have shown it is NP-Hard to 
test whether an instance of ACOPF is feasible. In our numerical experiments, 
however, we choose R such that each element is independently identically dis¬ 
tributed, uniformly over [0, 1], which need not be feasible for the instance of 
ACOPF. Subsequently, we compute t,h,g,u,v,z to match the R, projecting 
the values onto the intervals given by the box-constraints. Although one may 
improve upon this simplistic choice by a variety of heuristics, it still performs 
well in practice. 

The Choice of gk^Vk 

The choice of g,k,Vk may affect the performance of the algorithm. In order 
to prove convergence, one requires limfe_j.oo —>■ 0, but computationally, we 
consider 72. = 0, which obliterates the need for varying Vk and computation of 
the gradient, Vdet(i?^i?), as suggested by [T^ . 
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In order to prove convergence, any choice of > 0 is sufficient. Computa¬ 
tionally, we use /J,fe = /r = 0.0001 throughout much of the reported experiments. 
On certain well-known pathological instances, c.f.. Section 17.31 other choices 
may be sensible. Considering that many of those instances are very small, 
the use of large step-sizes, such as = 0.01 may seem justified. We 

have also experimented with an adaptive strategy for changing /x*,, as detailed 
in Section 17.41 There, we require the infeasibility to shrink by a fixed factor 
between consecutive iterations. If such an infeasibility shrinkage is not feasi¬ 
ble, one ignores the proposed iteration update and picks pfe+i < /ifc such that 
the infeasibility shrinkage is observed. Although the performance does improve 
slightly, we have limited the use of this strategy to Section mi We have also 
experimented with decreasing /j.^ geometrically, i.e. /Xfe = c • but we have 

observed no additional benefits. The fact one does not have to painstakingly 
tune the parameters is certainly a relief. 


5 Implementation Details 


In order to understand the details of the implementation, it is important to 
realise that 


• the quadratic forms such as Ykxx^ in the polynomial optimisation problem 
([PP4] I and Y^RR^ in the augmented Lagrangian (1551) are used for the 


simplicity of presentation. Indeed, the evaluation of the quadratic forms 
can be simplified considerably, when one does not strive for the brevity of 
expression 

many well-known benchmarks, including the Polish network, consider an 


extension of the polynomial optimisation problem ([PP4 


We elaborate upon these points in turn. 

5.1 The Simplifications 

Let us consider the example of triYkXx^) in more detail. The evaluation of the 
quadratic form can be simplified to 16 ||Tfe||o multiplications, i.e., performed in 
time linear in the number HlfeHo of non-zero elements of the matrix Yfc. More¬ 
over, the evaluation of the trace of the quadratic form requires only 2 HlfeHo 
multiplications, as one needs to consider only the diagonal elements. The num¬ 
ber of non-zero elements varies in each power system, and is quadratic in the 
number of buses for hypothetical systems, where there would be a branch be¬ 
tween each pair of buses, but the number of non-zero elements is linear in the 
number of buses, in practice. Consequently, one can simplify the first step of the 
inner loop (Lines 6-7), where one minimises the augmented Lagrangian (1331) . 
i.e., £{R,t,h,g,u,v,z,X^,X^,X^,X'^,X",X^), with respect to t, as follows. At 
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iteration k, for each coordinate j, in parallel, one computes the update: 


Ak) _ 


TT* + c]Sb + 2c]Sl{P^f' - 


A* 


i + HSI 


(35) 


where Sb is the base power of the per-unit system, such as 100 MVA, and tt* 
are the residuals: 

IT A .— . (oO) 

n 

Notice that the terms c]S'& and ‘2CjS^ are constant throughout the run and can 
hence be pre-computed, while the residuals (1361) have to be precomputed only 
once per iteration, e.g., just before the termination criteria are evaluated (Line 
10), and hence there are only 3 multiplications and 1 division involved, in ad¬ 
dition to the run-time of the evaluation of the quadratic form. Subsequently, 
one projects onto the box-constraints, i.e., if tj is large than the upper bound, 
it is set to the upper bound. If tj is smaller than the lower bound, it is set to 
the lower bound. Similar simplifications can be made for minimisation of the 
augmented Lagrangian (l33l) with respect to other variables. 


5.2 The Extensions 

In an often considered, but rarely spelled out [33] extension of the ACOPF 
problem (fT^ . one allows for tap-changing and phase-shifting transformers, as 
well as parallel lines and multiple generators connected to one bus. Let us 
denote the the total line charging susceptance (p.u.) by bi^, the transformer 
off nominal turns ratio by tim-, and the transformer phase shift angle by (pim- 
Then, the thermal limits (HU become: 


( cmax\2 

W/m ) 

-iv{ZimXx'^) 

-tv{ZlmXx'^) 



-tr{ZlmXx'^) 

1 

0 

^ 0 

(37) 

-tT{ZimXx'^) 

0 

1 



( cmax\2 

WZm ) 

-tr(T;^a;a;'^) 

-tr(T/„a;a;^) 



tr(Tj^a;x^) 

1 

0 

^ 0 

(38) 

tr(T/^a;x^) 

0 

1 




II 









where: 


ry 9im / r , 1 \ 

—-p^ye-iei + e;+|Ar|e;_,_|^| j 


(39) 


^Im 


+ 


9im cos(^/77T,) H- biYn COs(0/j77, H“ "2 ) 
gim H“ bim H" " 2 ) 


^ (e;e^ + e^ef + e/+|Ar|e^_,_|^| + em+\N\e-f^\N\) 


2 tlr. 


'^Im — 3tm(SmSrn + ^m+\N\£‘m+\N\) 

9lm COs( 0/m) “t“ ^/m COs( (j^lm "2 ) 


^ ' (eie^+|Ar| + e.m+\N\eJ — e-i+\N\e^ — e-meJj^\N\) 

(40) 

T : „ „T , „ T 


‘ 2 tlr, 


{eie^ + ernGj + ei+|Ar|e„_,_|^| + e^+|Ar|e,_,_|^|) 


9 lm sill( ^Im sill( (plm “t“ 2 ) 




2 bl m + bi r 


I ^LTTi / T I T \ 

— -- [eiBi + ei+|7v|e;+|jv|j 


2 t 


^ ' iBi+\N\em + - eje^_,_|^| - em+|jv|ef) 

(41) 


Im 


Qlm COs(^/77i) -|- biYTi COs((^;77T, -\- 2 ) 


+ 


'rim = - 


Qlm “h bi^fji ^) 

^bim H“ bi m ^ _ „T I „ „T \ 

^ K^m^rn ' ^m+\N\^m+\N\) 

9lm C0S( 4>lm) blm COs( (plm “I” "2 ) 
9lm sill( ^/m) H“ bi<fji sill( (plm "2 ) 


^ ^rn-\-\N\^l ^l+\N\^m 

(e/e^ + e-mef + e/+|Ar|e^+|Ar| + ^m-\-\N\G-'[^\N\) 


(42) 


^ (ez+|Ar|e^ + eme^|jv| - e/e^_,_|jv| - em+\N\ef) 
(eje^ + Smef + ej+|Ar|e^_,_|^| + em+\N\Bf+\N\) 


where Ck be the standard basis vector in K”, one can perform similar 
simplifications. 

For instance, following the pre-computation of vectors c^,c^,c^ € 
prior to the outer loop of the algorithm, the evaluation of the trace of the 
quadratic form, tr(ZimXx'^), considering (l39l) can be implemented using 4 look¬ 
ups into the vector x and 13 float-float multiplications: 

tr(Zi^xx^) = c^(xf + (43) 

^Im (2xiXm + 2a:;_|_|7v|2:m-|-|Af|) 

4” Blm.^^l^'m+\N\ 2XmXlJ^\J<[^. 


One can simplify the multiplication in a similar fashion for the remaining ex¬ 
pressions involving (traces of) quadratic forms of and T/m as well. 
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6 An Analysis 


The analysis needs to distinguish between optima of the semidefinite program¬ 
ming problem [LL] which is convex, and stationary points, local optima, and 
global optima of [RrBC] which is the non-convex rank-constrained problem. 


Let us illustrate this difference with an example: 

Example 6.1. Consider the following simple rank-constrained problem: 


max tr(diag(3,l)lT), (44) 

WeSl,rank(W) = l 


subject to ti{IW) = 1 and 

o) =(1,0)(1,0)^, 

^■=(o J) =(0,1)(0,1)^=:M^. 


W* is the unique optimal solution of (H41) , as well as the optimal solution of the 
SDP relaxation, where the rank constraint is dropped. There exists a Lagrange 
multiplier such that W is a stationary point, though. Let us define a Lagrange 
function £(i?. A) = tr(diag(3, l)i?i?'^)-|-A(tr(/i?i?^) —1). Then VrC{RR^, X) = 
2diag(3, 1)R + 2XR. If we plug in R and A = —1 we obtain VrC{RR^ , A) = 
2 diag(3,1)(0,1)^ -I- 2A(0,1)^ = (0,0)^. Hence (R, A) is a stationary point of 
a Lagrange function. However, one can follow the proof of Proposition 16.41 to 
show that [i?, 0] is not a local optimum solution of SDP relaxation. □ 


Indeed, in generic non-convex quadratic programming, the test whether a 
stationary point solution is a local optimum [45] is NP-Hard. For ACOPF, 
there are a number of sufficient conditions known, e.g. |43j . Under strong 
assumptions, motivated by Observation O 3 , we provide necessary and sufficient 
conditions based on those of Grippo et al. [20] . We use O for Kronecker product 
and Ir for identity matrix in 


Proposition 6.2. Consider the SDP relaxation obtained from [RrBC] by drop¬ 
ping the rank constraint and writing it down as min tr{QW) such that tr{AiW) = 
bi for constraints i = 1,..., m in variable W P 0,W € . If there exists an 

optimum W* with rank r for the SDP relaxation, for any point R G RR^ 

is a global minimiser of [RrBC] if and only if there exists a X* € R™ such that: 


Q 




i? = 0 


q + E\*a ^ 0 

i=l 


R^{A,®lr)R = b. 


(45) 


V* = 1, . . . , TO. 
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Proof. The proof is based on Proposition 3 of Grippo et al. [10] , which requires 
the existence of rank r optimum and strong duality of the SDP relaxation. The 
existence of an optimum of [LL] with rank r is assumed. One can use Theorem 
im and Theorem 1 of m to show that a ball constraint is sufficient for strong 
duality in the SDP relaxation, c.f. Observation O 4 . □ 


Notice that the test for whether there exists an optimum W* with rank r 
for the SDP relaxation is suggested Proposition 5, i.e. by solving for rank r + 1. 
Next, let us consider the convergence: 

Proposition 6.3. There exists an instance-specific constant r', such that for 
every r > r', whenever Algorithmic with G = det(i?^i?) and 

limfe^oo 0 produces solution with \im.k^ao — >■ 0 and a local opti¬ 
mum is generated within the inner loop (5-11), is an optimal 

solution to \[LL]\ Moreover, r' depends on the number of constraints m in the 
optimisation problem and is 0 {y/rn). 

Proof. The proof follows from Theorem 3.3 of Purer and Monteiro m- One 
can rephrase Theorem 3.3 to show that if G jg a bounded sequence 

such that: 


Cl. lim/j;_^oQ Tk 0 

02 : limfe —^00 = 0 

C 3 : liminffc_>.oo > 0 for all bounded sequences G 

X r 


C 4 : rank(i?’^^)) < r for all k 


every accumulation point of is an optimal solution of [LL] Let us 

show that these four conditions are satisfied, in turn. Condition Ci , which effec¬ 
tively says that W = R^^'> {R^^'> )'^ should be feasible with respect to constraints 
(1131®, is affected by the termination criteria of the inner loop, albeit only 
approximately for a finite k and finite machine precision. Conditions C 2 and 
C 3 follow from our assumption that R^^'> is a local optimum. The satisfaction 
of Condition C 4 can be shown in two steps: First, there exists an r', such that 
for every feasible semidefinite programming problem with m constraints, there 
exists an optimal solution with a rank bounded from above by r'. This r' is 
0{y/rn). This follows from Theorem 1.2 of Barvinok |4], as explained by Pataki 
[50]. Second, the TZ — det(i?^i?) regularisation forces the lower-rank optimum 
to be chosen, should there be multiple optima with different ranks. This can 
be seen easily by contradiction. Finally, one can remove the requirement on the 
sequence to be bounded by the arguments of Purer and Monteiro [H], as per 
Theorem 5.4. □ 


Further, 
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Proposition 6.4. Consider an instance of ACOPF such that there exists an 
optimum solution W* of \[LL]\ with rank 1 and'ik : c| > 0. Let he a iterate 
produced by the Algorithm[J\ when run with r = 1 and moreover it is such that 
Tk = 0. Then if RE?" (where R = is a local optimum of [R2BC], then 

j^(fc)(j^(fc))T ^ global optimum solution of [RIBC] and [PP4] 


Proof We will prove this proposition by contradiction. For the sake of contra¬ 
diction, we assume that i? is a local optimum and is not a global optimum. 
Therefore, we know that objective function of [RIBC] for Wi = W* is smaller 
than for W 2 = and both VFi and W 2 are feasible. By the assump¬ 

tion on optimum solution W* of ] [LL][ we know that W* can be written as W* = 
ww"’- . Now, it is easy to observe that if we define R^ = [(1 — \/X)R‘'^\ -v/Aw*] 
then for any A S [0,1] this vector is feasible for [R2BC]. Moreover, for A = 0 we 
have R^ = R. Because the objective function F of [RrBC] is convex in W with 
c| > 0, we have that 

F{R^{R^f) < (1 - X)F{W2) + XF{Wi) < F(W2) 

and hence for all A G (0,1] we have that F{LC'{R^Y') < F{W 2 ), which is a 
contradiction with the assumption that i? is a local optimum. □ 

Overall, 

Remark 1. The first iteration of the outer loop of Algorithm [T] may produce a 
global optimum to |[Rl]| and [RrBC] and [PP4] as suggested in ProDOsitions l3.ll - 
16.41 The second and subsequent iterations of the outer loop of Algorithm [T] may 
find the optimum of |[LL][ the semidefinite programming problem, as suggested 
in Proposition 16.31 but for 72. = 0, they are not guaranteed to find the global 

^[PPI 


optimum of [R1 


nor 


RrBC] 


nor 


One should also like contrast the ability to extract low-rank solutions with 
other methods: 


Remark 2. Whenever there are two or more optima of [LL] with two or more 
distinct ranks, the maximum rank solutions are in the relative interior of the 
optimum face of the feasible set, as per Lemma 1.4 in [28]. Primal-dual interior- 
point algorithms for semidefinite programming, such as SeDuMi |55j . in such a 
case return a solution with maximum rank. 

Put bluntly, one may conclude that as long as one seeks the exact optimum 



of whether one uses 72 = det(7?^i?) or 0. 


Although 72 = 0 or z/ = 0 does 
in some cases lai. 


not guarantee the recovery of low-rank solutions of [LL] 

72 = det(7?^7?) and > 0 does not guarantee that the low-rank solution of the 
augmented Lagrangian (l33l) coincides with the optimum of [ACOPFH in some 
other cases. It may hence be preferable to use 72 = 0, as it allows for more 
iterations per second and post hoc testing of global optimality. 
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7 Numerical Experiments 

We have implemented Algorithm [1] in C++ with GSL and OpenMP and tested 
it on a collection |62j of well-known instances. For comparison, we have used 
three solvers specialised to ACOPF: 

• the MATLAB-based Matpower Interior Point Solver (MIPS) version 1.2 
(dated March 20th, 2015), which has been developed by Zimmerman et 

al. [02] 

• the C-based Semidefinite Programming optimal Power Flow (sdp_pf) ver¬ 
sion 1.0 (dated January 17th, 2014), which has been developed by Mohlzahn 
et al. [44] using SeDuMi of Sturm et al. [55] 

• the MEX-based 0PF_Solver beta version dated December 13th, 2014, 
listed as the most current as of June 1st, 2 oid!l, which has been developed 
by Lavaei et al. [35] using and SDPT3 of Tutiincii et al. [S01|S7]. Notice 
that 0PF_Solver produces feasible points with objective values that are 
near the values of the global optima across all instances tested, for some, 
non-default settings; we have used the per-instance settings of epB, epL, 
and line_prob, as suggested by the authors. 

For the comparison presented in Table [U we have used a standard laptop with 
Intel i5-2520M processor and 4 GB of RAM. We believe this is fair, as the solvers 
we compare with cannot make a good use of a more powerful machine. For the 
presentation of scalability of our code, we have used a machine with 24 Intel 
E5-2620v3 clocked at 2.40GHz and 128GB of RAM, but used only the numbers 
of cores listed. 

We have also tested six general-purpose semidefinite-programming (SDP) 
solvers: 

• CSDP version 6.0.1, which has been developed by Borchers [7] 

• MOSEK version 7.0, which has been developed by MOSEK ApS [Tl IST] 

• SeDuMi version 1.32, which has been developed by Sturm et al. [55] 

• SDPA version 7.0, which has been developed by the SDPA group m 

• SDPT3 version 4.0, which has been developed by Tutiincii et al. [5i[57] 

• SDPLR version 1.03 (beta), which has been developed by Burer and Mon- 
teiro [Tl] 

Five of the codes (GSDP, SeDuMi, SDPA, SDPLR, and SDPT3) have been 
tested at the NEOS 7 facility m at the University of Wisconsin in Madison, 
where there are two Intel Xeon E5-2698 processors clocked at 2.3GHz and 192 
GB of RAM per node. MOSEK has been used at a cluster equipped with one 
AMD Opteron 6128 processor clocked with 4 cores at 2.0 GHz and 32 GB of 
RAM per node, as per the license to Lehigh University. 

^http://ieor.berkeley.edu/-lavaei/Software.html 
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7.1 IEEE Test Cases 


Our main focus has been on the IEEE test cases. In Table [I] we compare the 
run-time of our implementation of Algorithm [T] with the run-time of the three 
leading solvers for the ACOPE listed above, two of which (sdp_pf, 0PF_Solver) 
use elaborate tree-width decompositions. In order to obtain the numbers, we 
ran Matpower first using default settings, record the accuracy with respect to 
squared error Tk (IMl) . ran our solver up to the same accuracy, and record the 
time and the objective function. 

In Table [21 we present the run-time of six popular general-purpose SDP 
solvers for comparison. We note that we have used the default parameters and 
tolerances of each code, so the precision may no longer match Matpower. We 
list the reported CPU time rounded to one decimal digit, if that yields a non¬ 
zero number, and to one significant digit, otherwise. When we display a dash, 
no feasible solution has been found; the severity of the constraint violation and 
abruptness of the termination of the solver vary widely. For example, MOSEK 
terminates very abruptly with fatal error stopenv on five of the instances, while 
SDPT3 often runs into numerical difficulties. Outside of SDPLR on the largest 
instance, no solver ran out of memory, iteration limit, or time limit. In this 
comparison, SeDuMi seems to be most robust solver. SDPLR is also rather 
robust, but three orders of magnitude slower. Throughout, all interior-point 
methods (CSDP, Mosek, SeDuMi, SDPA, and SDPT3) perform much better on 
the dual of the SDP, than on the primal. Please note that direct comparison 
with the run-time of our implementation of Algorithm |T] reported in Table [I] is 
no possible, due to the use of three different platforms: Five of the codes (CSDP, 
SeDuMi, SDPA, SDPLR, and SDPT3) have been tested at the NEOS facility, 
which does not allow for our code to be run, while Mosek has been tested at a 
Lehigh University facility. Still, either facility has machines considerably more 
powerful than the machine used in the tests above, and we report the CPU 
time reported by the individual solvers, rather than the wall-clock time, so we 
believe it is fair to claim our specialised solver outperforms the general-purpose 
solvers. 

7.2 NESTA Test Cases 

Next, we have tested our approach on the recently introduced test cases from 
the NICTA Energy System Test Case Archive (NESTA) [15]. There, bounds on 
commonly used IEEE test cases have been carefully tightened to make even the 
search for a feasible solution difficult for many solvers. For example, in the so 
called Active Power Increase (API) test cases, the active power demands have 
been increased proportionally throughout the network so as to make thermal 
limits active. In Table jS] we present the results. Out of the 11 API instances 
tested, Matpower and sdp_pf fail to find feasible solutions for three instances 
each. Our implementation of Algorithm [l] using default settings, including 
TZ = 0, obtains feasible solutions across all the 11 instances, while improving 
over the objective function values obtained by either Matpower or sdp_pf. For 
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example, on the instance nesta_case30_as__api, we improve the objective function 
value by about 10%, while on the instance nesta_case30_fsr__api, we improve 
the objective by considerably more than 10%. Although we have not tested 
all instances of the NESTA archive, we are very happy with these preliminary 
results. 


7.3 Pathological Instances 


Additionally, we have tested our approach on a number of recently introduced 
pathological instances [3Tl|45l|9l[T9l[25]. Depending on the choices of R, v, and 
Algorithm [T] may perform better or worse than the relaxation of Lavaei and Low 
([LL] I. Let us illustrate this on the suggested settings of 72. = 0, arbitrary, 
and single-thread execution. In the example of Bukhsh et al. [5], known as 
case2w, there are two local optima. Whereas many heuristics may fail or find 
the local optimum with cost 905.73, we are able to find the exact optimum with 
cost 877.78 in 0.49 seconds up to the infeasibility of 9.38 • 10“® with ^ = 0.01 
and up to infeasibility of 7.23 • 10“® in 0.94 seconds with /r = 0.0001. When 
we replace = 1.05 with 1.022, as in [19], the instance becomes harder 

still, and the optimum of the relaxation of Lavaei and Low ([LL] I is not rank- 
1 . (Although one could project onto the feasible set of [RIBC] and extract 
a feasible solution of [PP4][ it would not be the optimum of [PP4] ) With 


fi = 0.01, we find only a local optimum with cost 888.05 up to infeasibility 
of 8.29 • 10“®, but with /r = 0.0001, we do find the local optimum, which in 
our evaluation has cost 905.66 up to infeasibility of 1.46 • 10“®. We stress that 
this optimum is not the solution of plain ([LL]), but at the same time that. 


we do not provide any guarantees of improving upon the relaxation of Lavaei 
and Low ([LL] |. In the example case9mod of Bukhsh et al. [9], ^ = 0.01 does 


not converge within 10000 iterations, but using /r = 0.0001, we are able to hnd 
the exact optimum of 3087.84 and infeasibility 9.43 • 10“^^ in 12.16 seconds. 
In the example case39mod2 of Bukhsh et al. [S], we are able to find only a 
local optimum with cost 944.71 and infeasibility 5.63-10“® after 45.03 seconds, 
whereas the present-best known solution has cost 941.74. In the example of 
Molzahn et al. EH ESI, which is known as LMBM3, we find a solution with 
cost 5688.10 up to infeasibility of 9.98 • 10“® in 0.88 seconds using /r = 0.01 
and a solution with cost 5694.34 up to infeasibility of 1.14 • 10“® in 0.56 using 
/r = 0.0001. This illustrates that the choice of /r is important and the default 
value, albeit suitable for many instances, is not the best one, universally. 
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Table 1: The results of our numerical experiments on IEEE test cases and certain well-known pathological instances. Dash (-) 
indicates no feasible solution has been provided, often due to numerical issues or the lack of convergence. 


Instance 

MATPOWER 

sdp pf 

0PF Solver 

Algorithm [T] 

Name 

Ref. 

Obj. 

Time [s] 

Obj. 

Time [s] 

Obj. 

Time [s] 

Obj. 

Time [s] 

case2w 

El 

— 

— 

877.78 

17.08 

877.78 

2.52 

877.78 

0.077 

caseSw 

El 

— 

— 

560.53 

0.56 

560.53 

2.70 

560.53 

0.166 

case5 

m 

1.755e-f04 

21.80 

1.482e-f03 

120.73 

1.482e-F03 

25.93 

1.482e-F03 

0.263 

casebww 

[59] 

3.144-F03 

0.114 

3.144e-f03 

0.74 

3.144e-F03 

2.939 

3.144e-F03 

0.260 

casel4 

m 

8.0826-|-03 

0.201 

8.082g-1-03 

0.84 

- 

- 

8.082e-F03 

0.031 

caseSO 

m 

5.769e-f02 

0.788 

5.769e-f02 

2.70 

5.765e-F02 

6.928 

5.769e-F02 

0.074 

case39 

m 

4.189e-f04 

0.399 

4.189e-f04 

3.26 

4.202e-F04 

7.004 

4.186e-F04 

0.885 

case57 

m 

4.174e-f04 

0.674 

4.174e-f04 

2.69 

- 

- 

4.174e-F04 

0.857 

case57Tree 

m 

12100.86 

6.13 

* 1.045e-f04 

4.48 

* 1.046e-F04 

342.00 

* 1.046e-F04 

3.924 

easel 18 

m 

1.297e-f05 

1.665 

1.297e-f05 

6.57 

- 

- 

1.297e-F05 

1.967 

caseSOO 

m 

7.197e-f05 

2.410 

- 

17.68 

- 

- 

7.197e-F05 

90.103 


*: We note that the precision here is approximately 10 which is the case for all three SDP-based solvers, Algorithm [U as 
well as sdp_pf and DPF_Solver. 













Table 2: For comparison, the CPU time in seconds of certain well-known general-purpose SDP solvers on the Lavaei-Low 
primal or dual SDP obtained from IEEE test cases and certain well-known pathological instances. Dash (-) indicates no 
feasible solution has been provided, often due to numerical issues or the lack of convergence. 


Instance 



Time 

[s] 



Name 

SDP 

SeDuMi 1.32 

SDPA 7.0 

SDPLR 1.03 

CSDP 6.0.1 

SDPT3 4.0 

Mosek 7.0 

case2w 

primal 

0.3 

0.003 

0 

0.02 

0.4 

0.1 

case2w 

dual 

0.3 

0.004 

0 

0.02 

0.4 

0.1 

case5w 

primal 

0.4 

- 

0 

0.04 

0.5 

0.2 

case5w 

dual 

0.4 

- 

0 

0.03 

0.5 

0.2 

case9 

primal 

1.0 

0.05 

25 

0.1 

0.5 

0.3 

case9 

dual 

1.0 

0.07 

25 

0.1 

0.6 

0.3 

caseI4 

primal 

0.7 

0.05 

58 

0.1 

- 

0.2 

caseI4 

dual 

0.7 

0.04 

17 

0.1 

- 

0.3 

caseSO 

primal 

2.8 

- 

829 

0.8 

- 

- 

caseSO 

dual 

6.1 

- 

282 

1.0 

- 

- 

caseSOQ 

primal 

2.8 

- 

831 

0.8 

- 

- 

caseSOQ 

dual 

6.2 

- 

195 

1.0 

- 

- 

case39 

primal 

4.4 

- 

2769 

1.0 

- 

0.7 

case39 

dual 

7.7 

- 

723 

1.3 

- 

- 

case57 

primal 

3.2 

- 

1930 

0.5 

- 

0.7 

case57 

dual 

4.0 

- 

1175 

1.0 

- 

1.8 

easel18 

primal 

10.3 

0.9 

4400 

3.4 

- 

1.7 

easel18 

dual 

17.1 

- 

- 

13.0 

- 

- 

case300 

primal 

27.5 

1.7 

- 

109.6 

- 

- 

case300 

dual 

133.7 

- 

- 

66.7 

- 

- 



1600 



Figure 1: The evolution of the objective function value over time on instance 
case5, for a number of choices of and the adaptive updates of /i. 

7.4 Adaptive Updates of // 

As it has been stressed above, it is important to pick an appropriate fi. Alter¬ 
natively, one can adapt a strategy to change adaptively. This can be achieved 
by requiring the infeasibility to shrink by a fixed factor between consecutive 
iterations. If such an infeasibility shrinkage is not feasible, one ignores the pro¬ 
posed iteration update and decreases fi. Figures [1] and O show the evolution of 
objective function value and the infeasibility for 400 iterations with fixed values 
/r £ {10“^, 10“^, 10“^, 10“"^, 10“®, 10“®} and for the adaptive strategy on the 
instance case5. Figure [3] details the evolution of ^ within the adaptive strategy 
on the same instance. One can see that choosing very small /r forces the algo¬ 
rithm to “jump” close to a feasible point, and get “stuck”. On the other hand, 
larger fixed values of ^ lead to the same objective value. One can also see that 
the adaptive strategy is rejecting the first 80 updates, until the value of fi is 
small enough; subsequently, it starts to make rapid progress. 

7.5 Large Instances 

We have also experimented with the so called Polish network [62], also known 
as case2224 and case2736sp, as well as the instances collected by the Pegase 
project [SS], such as case9241peg. There, we cannot obtain solutions with the 
same precision as Matpower, within comparable run-times. However, we can 
decrease the initial infeasibility by factor of 10^ for case2224 in 200 seconds, by 
factor of 10® for case2736sp in 200 seconds, and by factor of 10"^ for case9241peg. 
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Table 3: The results of our numerical experiments on NESTA test cases. When no feasible solution has been provided, we list 
the infeasibility of the least infeasible point provided in red parentheses. 


Instance 

MATPOWER 

sdp pf 

Algorithm [T] 

Name 

Obj. 

Time [s] 

Obj. 

Time [s] 

Obj. 

Time [s] 

nesta_case3_lmbd 

5.812e+03 

0.946 

5.789e+03 

2.254 

5.757e+03 

0.149 

nesta_case4_gs 

1.564e+02 

1.019 

1.564e+02 

2.392 

1.564e+02 

0.139 

nesta_case5_pjm 

(1.599e-01) 

0.811 

(1.599e-01) 

2.708 

2.008e+04 

0.216 

nesta_case6_c 

2.320e+01 

0.825 

2.320e+01 

2.392 

2.320e+01 

0.379 

nesta_case6_ww 

3.143e+03 

0.884 

3.143e+03 

2.776 

3.148e+03 

0.242 

nesta_case9_wscc 

5.296e+03 

1.077 

5.296e+03 

2.621 

5.296e+03 

0.211 

nesta_casel4_ieee 

2.440e+02 

0.762 

2.440e+02 

3.164 

2.440e+02 

0.267 

nesta_case30_as 

8.031e+02 

0.861 

8.031e+02 

3.916 

8.031e+02 

0.608 

nesta_case30_fsr 

5.757e+02 

0.898 

5.757e+02 

6.201 

5.750e+02 

0.613 

nesta_case30Jeee 

2.049e+02 

0.818 

2.049e+02 

4.335 

2.049e+02 

0.788 

nesta_case39_epri 

9.651e+04 

0.863 

9.649e+04 

5.676 

9.651e+04 

0.830 

nesta_case57 Jeee 

1.143e+03 

1.119 

1.143e+03 

8.053 

1.143e+03 

0.957 

nesta_casell8Jeee 

(4.433e+02) 

1.181 

(4.433e+02) 

18.097 

4.098e+03 

4.032 

nest a.case 162 Jeee.dtc 

(3.123e+03) 

0.975 

(3.123e+03) 

32.153 

4.215e+03 

8.328 

nesta_case3_lmbd__api 

3.677e+02 

1.113 

3.333e+02 

2.459 

3.333e+02 

0.029 

nesta_case5_pjm__api 

(5.953e+00) 

1.002 

(5.953e+00) 

2.426 

3.343e+03 

1.059 

nesta_case6_c__api 

8.144e+02 

1.034 

8.144e+02 

2.456 

8.139e+02 

0.722 

nesta_case9_wscc__api 

6.565e+02 

1.135 

6.565e+02 

2.869 

6.565e+02 

0.708 

nesta_casel4Jeee__api 

3.255e+02 

1.035 

3.255e+02 

3.569 

3.245e+02 

0.214 

nesta_case30_as__api 

5.711e+02 

1.156 

5.711e+02 

5.433 

5.683e+02 

0.615 

nesta_case30_fsr__api 

3.721e+02 

0.862 

3.309e+02 

5.276 

2.194e+02 

0.604 

nesta_case30Jeee__api 

4.155e+02 

1.076 

4.155e+02 

5.293 

4.155e+02 

0.618 

nesta_case39_epri__api 

(1.540e+02) 

0.965 

(1.540e+02) 

5.109 

7.427e+03 

1.978 

nesta_case57 Jeee__api 

1.430e+03 

1.041 

1.429e+03 

7.625 

1.429e+03 

2.437 

nest a.case 162 Jeee_dtc__api 

(1.502e+03) 

1.215 

(1.502e+03) 

43.339 

6.003e+03 

5.833 












Figure 2: The evolution of the infeasibility over time on instance case5, for a 
number of choices of /r and the adaptive updates of /i. 



Figure 3: The evolution of /i with adaptive updates over time on instance caseS. 
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Figure 4: The speed-up as a function of the number of cores employed. 


again in 200 seconds. Although these results are not satisfactory, yet, they may 
be useful, when a good initial solution is available. As discussed in Section El 
we aim to improve upon these results by using a two-step method, where first- 
order methods on the convex problem are combined with a second-order method 
on the non-convex problem. These results also motivate the need for parallel 
computing. 

7.6 Parallel Computing 

As it has been discussed above. Algorithm [T] is easy to implement in a parallel 
fashion. We have implemented our multi-threading variant using OpenMP, in 
order to achieve portability. The significant overhead of using OpenMP makes it 
impossible to obtain a speed-up on very small instance. For example, in caseSw, 
there is not enough work to be shared across (any number of) threads to offset 
the overhead. In Figure IH we present the scaling properties of Algorithm [TJ 
In particular, we illustrate the speed-up (the inverse of the multiple of single- 
threaded run-time) as a function of the number of cores used for a number of 
larger instances. We can also observe that even easel 18 is too small to benefit 
from a considerable speed-up. On the other hand, for large datasets, such as 
the Polish case2224 or larger, we observe a significant speed-up. 
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8 Conclusions 


Our approach seems to bring the use of SDP relaxations of ACOPF closer to the 
engineering practice. A number of authors have recently explored elaborately 
constructed linear programming and second-order cone programming Emu 
SSIIMIIIS] relaxations, many EZIES] of which are not comparable to the SDP 
relaxations in the sense that they may be stronger or weaker, depending on the 
instance, but aiming to be solvable faster. Algorithm [T] suggests that there are 
simple first-order algorithms, which can solve SDP relaxations of ACOPF faster 
than previously, at least on some instances. 

A major advantage of first-order methods over second-order methods [621155] 
is the ease of their parallelisability. This paper presents considers a symmetric 
multi-processing, where memory is shared, but a distributed variant, where C 
agents perform the iterations, is clearly possible, c.f. [41]. The C agents may 
represent companies, each of whom owns some of the generators and does not 
want to expose the details, such as cost functions. If C agents are C computers, 
a considerable speed-up can be obtained. Either way, power systems analysis 
could benefit from parallel and distributed computing. 

The question as to how far could the method scale, remains open. In order 
to improve the scalability, one could try to combine first- and second-order 
methods [33]. We have experimented with an extension, which uses Smale’s 
a-P theory [6] to stop the computation at the point zq, where we know, based 
on the analysis of the Lagrangian and its derivatives, that a Newton method 
or a similar algorithm with quadratic speed of convergence m will generate 
sequence Zi to the correct optimum z*, i.e. 

\z,-z*\<il/2f-^\zo-z*\. (46) 

This should be seen as convergence-preserving means of auto-tuning of the 
switch to a second-order method for convex functions. 

One should also study infeasibility detection. Considering the test of feasi¬ 
bility of an SDP is not known to be in NP [53], we assumed the instance are 
feasible, throughout. Nevertheless, one may need to test their feasibility first, in 
many practical applications. There are some very fast heuristics [47] available 
already, but one could also use Lagrangian methods in a two-phase scheme, 
common in robust linear-programming solvers. Both in methods based on sim¬ 
plex and feasible interior-point, one first considers a variant of the problem, 
which has constraints relaxed by the addition of slack variables, with objective 
minimising a norm of the slack variables. This makes it possible to find feasible 
solutions quickly, and by duality, one can detect infeasibility. 

One could also apply additional regularisations, following [MUSS]. In |[Rl][ 
one could drop the rank-one constraint and modify the objective function to 
penalise the solutions with large ranks, e.g., by adding the term A||VF||* to the 
objective function, where || • ||* is a nuclear norm [18] and A > 0 is a parameter. 
Alternatively, one can replace the rank constraint by the requirement that the 
nuclear norm of the matrix should be small, i.e. ||IF||, < A. However, both 
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approaches require a search for a suitable parameter A such that the optimal 
solution has indeed rank 1. Moreover, the penalised alternative may not produce 
an optimal solution of [PP4] necessitating further algorithmic work. 


Finally, one could extend the method to solve relaxations a variety of related 
applications, such as security-constrained problems, stability-constrained prob¬ 
lems, network expansion planning [40) . and unit commitment problems. The 
question as to whether the method could generalise to the higher-order relax¬ 
ations, c.f., Ghaddar et al. m, also remains open. First steps m have been 
taken, but much work remains to be done. 
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